Image charge effects in single-molecule junctions: Breaking of symmetries and 
negative differential resistance in a benzene transistor 
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Both experiments and theoretical studies have demonstrated that the interaction between the cur- 
rent carrying electrons and the induced polarization charge in single-molecule junctions leads to a 
strong renormalization of molecular charging energies. However, the effect on electronic excitations 
and molecular symmetries remain unclear. Using a theoretical framework developed for semiconduc- 
tor nanostructure based single-electron transistors (SETs), we demonstrate that the image charge 
interaction breaks the molecular symmetries in a benzene based single-molecule transistor operating 
in the Coulomb blockade regime. This results in the appearance of a so-called blocking state, which 
gives rise to negative differential resistance (NDR). We show that the appearance of NDR and its 
magnitude in the symmetry-broken benzene SET depends in a complicated way on the interplay 
between the many-body matrix elements, the lead tunnel coupling asymmetry, and the bias polar- 
ity. In particular, the current reducing property of the blocking state causing the NDR, is shown to 
vanish under strongly asymmetric tunnel couplings, when the molecule is coupled stronger to the 
drain electrode. The calculated IV characteristic may serve as an indicator for image charge broken 
molecular symmetries in experimental situations. 
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I. INTRODUCTION 

Over the past decade the field of single-molecule elec- 
tronics has experienced significant progress both ex- 
perimentally and theoretically. On the experimen- 
tal side break-junction and electromigration techniques 
have become standard ways of realizing nano-scale junc- 
tions with only one molecule bridging the gap be- 
tween the source and drain electrodes With the cou- 
pling between the molecule and the electrodes often 
being a highly uncontrollable parameter, the transport 
through single-molecule junctions spans different param- 
eter regimes depending on the fabrication techniques. 4 
Three-terminal measurements where the molecule cou- 
ples capacitively to a third gate electrode most often end 
up in the weak coupling regime where the Coulomb block- 
aded transport through the molecule is dominated by se- 
quential tunneling. In these setups, the coupling to the 
gate electrode allows control of the alignment between 
the molecular levels and the Fermi levels of the leads. 
The IV characteristics is therefore often summarized in 
the so-called charge stability diagram, which maps out 
the differential conductance dI/dV s d as a function of the 
gate and source-drain voltage. These diagrams help to 
understand the excitations intrinsic to the molecule, as 
for example the vibrational degrees of freedomj2£ 

One still unresolved issue is how the junction polar- 
ization resulting from the charging of the molecule in- 
fluences the molecular states. In, for example, three- 
terminal electromigrated nanoscale junctions,—^ where 
the dimensions of the source and drain electrodes are 
large compared to the metallic screening length, this 



effect can be expected to be significant. So far, 
both experiments^ and theoretical simulations^ have 
demonstrated that the interaction with the polarization 
charge — the so-called image charges — leads to a strong 
renormalization of the molecular charging energies. In a 
charge-stability diagram, this is reflected in addition en- 
ergies reduced with up to several electron-volts compared 
to the value expected from gas phase level positions of the 
molecule. On the other hand, the effect on excited states 
remain unresolved. A number of experimental studies 
have indicated image charge stabilized states close to the 
electrodes and modified excitation energies.^ 1 ^ Hence, a 
thorough investigation of the image charge effect and its 
influence on model parameters can help to improve the 
interpretation of experimental observations. 

With the recent interest in the role of molecular de- 
generacies in the Coulomb blockade regime* 1 -^ 13 it is 
of relevance to address the effect on such degeneracies. 
While splittings of degenerate levels on the order of the 
tunnel broadening results in quasi-degenerate states that 
still behave as being degenerate, a larger splitting com- 
pletely destroys interference effects and the associated 
signature in the transport characteristics . 11 ' 13 Degenera- 
cies also play an important role in Jahn- Teller active 
molecules that undergo distortions upon charging. Here, 
the higher-dimensional adiabatic potential energy surface 
of the charged molecule resulting from the coupling to the 
Jahn- Teller active vibrations leads to distinct transport 
characteristics 1 14 > 15 Also in this case can a splitting of the 
degenerate levels result in a qualitatively different signa- 
ture in the stability diagram. For example, a transport 
signature characteristic of the pseudo Jahn- Teller effect 
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might result when the size of the energy splitting matches 
a multiple of a vibrational energy J£ 

In this paper, we present a quantitative study of the 
impact on molecular degeneracies by studying the image 
charge effect in a benzene-based single-molecule junction. 
Despite the presence of the Jahn- Teller effect in the ben- 
zene moleculejil we will focus on the splitting of the elec- 
tronic states due to image charge effects and neglect any 
coupling to vibrations. 

Recent theoretical studies have already touched upon 
the image charge effect and other symmetry-breaking 
agents in a benzene SET, however, without taking into 
account the full interaction between the molecule and the 
image charges! 11 ' 18 Here, the purpose is twofold. First, 
we present a general framework that has previously been 
used to account for image charge effects in semiconductor 
nanostructures^ and discuss its applicability for single- 
molecule junctions operating in the Coulomb blockade 
regime. Second, the implication on the degeneracies 
of the benzene molecule and the consequence for the 
low-bias transport characteristics are investigated. We 
find that the image charge effect indeed breaks the high 
symmetry of the molecule and leads to a large splitting 
(~ 40 — 80 meV) of the degenerate ground-state of the 
singly-charged molecule. 

Due to a breakdown of the transport selection rules 
that apply in the isolated molecule, the symmetry-split 
excited state of the charged molecule is turned into a 
so-called blocking state. As a consequence, an NDR fea- 
ture appears at a bias corresponding to the level split- 
ting. The stability of the NDR feature with respect to 
an asymmetry in the lead couplings and the coupling 
site on the molecule is analyzed. Similar NDR features 
caused by radiative relaxation to a blocking stated and 
an interference-induced blocking stat e 11 ' 21 have been re- 
ported previously in the literature for the benzene SET. 
Contrary to these cases, the NDR feature observed here 
is exclusively caused by the broken symmetry of the 
molecule and needs none of these additional effects to 
occur. 

The paper is organized as follows. Sec. [TTI presents the 
theoretical approach used to describe the image charge 
effect in this work. This involves a generalized Hamilto- 
nian for the molecule with extra terms originating from 
the interaction with the image charges. In Sec. IIII1 the 
benzene SET and the semi-empirical Pariser-Parr-Poplc 
Hamiltonian is introduced. Furthermore, an analysis of 
the broken symmetry and the accompanying splitting of 
the degenerate molecular states together with a condi- 
tion for the occurrence of NDR is given. The resulting 
low-bias IV characteristics are presented and analyzed 
in Sec. IIV1 Finally, Sec. IVl summarizes our findings and 
points to other situations where the image charge effect 
may be important in single-molecule junctions. 



II. JUNCTION HAMILTONIAN AND THE 
CURRENT 

In the following section, we present a general frame- 
work for the description of single-molecule junctions op- 
erating in the Coulomb blockade regime, i.e. with a 
weak tunnel coupling between the molecule and the leads. 
While different approaches to the calculation of the cur- 
rent to different orders in the leads couplings have been 
given elsewhere,— the present study focuses on the in- 
teraction with the image charges in the junction environ- 
ment. Figure QJ, illustrates an idealized setup for such a 
three-terminal junction with the molecule lying on a gate 
oxide between the source and drain electrodes. With the 
molecule playing the role of the quantum dot, this type of 
single- molecule junctions is very similar to multi-terminal 
semiconductor-based quantum dot structures, where the 
theoretical description of image charges and their influ- 
ence on the dot states is well-established^ The theo- 
retical framework is based on a macroscopic description 
of the metallic electrodes and the surrounding dielectric, 
while the quantum dot(s) is described quantum mechan- 
ically. 

It is here illustrated how the same formalism can be ap- 
plied to single-molecule junction. We briefly summarize 
the most important aspects of the approach and explain 
in detail the appearance of new terms in the molecu- 
lar Hamiltonian that result from the image charge effect. 
The validity of a macroscopic electrostatic description for 
single-molecule SETs is also discussed. 

A. Hamiltonian 

A general Hamiltonian for a single-molecule junction 
including the image charge effect, can be written as the 
sum of the following terms, 

H = H mo \ + Ht + H i oa ds + ifmol-cnv + H cnv . (1) 

Here, ff mo i denotes the Hamiltonian of the molecule ex- 
pressed in a basis of atomic orbitals 4>i with correspond- 
ing creation and annihilation operators {c\ a , c ir7 }. In the 
present work, this basis will be the atomic p z orbitals 
of the carbon atoms in the benzene molecule. The next 
two terms account for the tunnel couplings to the lead 
(labeled by a) electrons, 

H T + Pleads = C4<7« C aa + t-C + ^ £ fc<r c La C ka a ■ 

kaa kaa 

(2) 

For simplicity, it is here assumed that only one orbital 
<j> a on the molecule couples to each of the left (a = L) 
and right (a — R) leads. 

The remaining two terms in Eq. ([T]) account for the 
electrostatic energy arising from the interaction between 
the charge of the molecule and junction. This includes 
the interaction with the induced polarization charge in 
the environment due to charging of the molecule (the 
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a) 



b) 



• para 

• meta 



c) 



d) 




FIG. 1: (Color online) Schematic illustration of a benzene single-electron transistor, a) Side view of the nanojunction with the 
molecule positioned on the gate dielectric between the source and drain electrodes. The capacitive coupling between the gate 
electrode (gray) at the bottom and the molecule, provides control over the energy levels of the molecule, b), c) and d) Top 
view of the junction with the molecule positioned in different configuration. The red and green dots indicate the coupling sites 
for coupling in the para and meta configuration, respectively, b) Symmetric setup with the molecule positioned in the middle 
of the gap between the electrodes, c) Asymmetric setup with the molecule positioned closer to the left electrode, d) Same as in 
c) but with an additional rotation of the molecule. Due to the interaction with the image charges of the junction, the rotation 
of the molecule breaks all the symmetries of the molecular Hamiltonian and results in qualitatively different IV-characteristics 
and NDR. 



image charges), the energy of the molecule in the poten- 
tial from voltages V* applied to the electrodes, and the 
electrostatic energy of the capacitively coupled electrodes 



J 



of the junction. For now, these two terms are most con- 
veniently expressed together with the molecular Hamil- 
tonian -ff mo i in real-space and first-quantization as 19 



#mol + H, 



mol-' 



1 N 
2m ^ 



N 
n=l 



V 



N 



tEE G(r„,rv)+5>C^, 



(3) 



n— 1 n'y^n 



where m is the electron mass, N is the number of elec- 
trons N on the molecule, G denotes the electrostatic 
Green's function of the junction, CV, is the capacitance 
matrix of the junction coupling the charge on the differ- 
ent electrodes. The single-particle potential V is given 
by 



V(t) 



dr'G(r,r')p ion (r')+J2<Xi(r)Vi, 



( 4 ) 

where G is the part of the Green's function that accounts 
for the induced potential due to an electron in r (see be- 
low), pi on is the charge distribution of ionic cores of the 
molecule (often described by a pseudo-potential in ab- 
initio approaches), and on is an electrode specific func- 
tion, which gives the spatial profile of the junction po- 
tential with a unit voltage applied to the i'th electrode 
and in the absence of the molecule. The cv-functions are 
solutions to Laplace equation 

- V- [e £ r (r)V a ,(r)] =0 (5) 

with boundary conditions Vj — Sij at the electrodes, e r is 



the spatially varying dielectric constant of the junction, 
and eo is the vacuum permittivity. 

Apart from the capacitive energy in last term of Eq. §3§ 
and the contribution from the applied voltages in the 
last term of Eq. the Hamiltonian in Eq. ^ bears 
close resemblance to the usual many-body Hamiltonian 
for interacting electrons moving in an external potential. 
However, the electrostatic Green's function plays an im- 
portant role in this modified Hamiltonian. It is seen to 
have replaced the Coulomb interaction (third and second 
terms of Eqs. (J3|) and respectively) and results in an 
additional single-particle term (first term of Eq. (Q|). 

The electrostatic Green's function solves Poisson's 
equation with a (^-function source term, 



V-[e er(r)VG(r,r')] = <J(r - r' 



(6) 



with Dirichlet boundary conditions Vj = 0. For a given 
junction geometry the Green's function gives the poten- 
tial in r due to a unit charge in r'. When both r and r' 
belong to the vacuum region of the junction where the 
molecule resides, the Green's function can be written as a 
sum of the direct (unscreened) Coulomb interaction plus 
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a contribution G from the image charges (screening), 

G{v,v') = y^- + G{v,v'). (7) 

The replacement of the Coulomb interaction with the 
electrostatic Green's function in Eqs. ([3]) and ((H) there- 
fore corresponds to a screening of the Coulomb inter- 
action by the polarization response of the junction (the 
image charges). The additional single-particle term in 
Eq. (0| has the form of an electronic self-interaction given 
by the G-part of the Green's function. This is the energy 
of the electron in its own induced potential, i.e. the in- 
teraction between the electron and its own image charge. 
For the standard example of a point charge positioned 
at a distance z from an infinite conduction surface this 
term reduces to the classical — l/4z-energy well known 
from classical electrostatics^ 

In the present work, the electrostatic Green's func- 
tion is obtained for the simplified junction geometry dis- 
cussed in App. [X] This allows for an analytical solution 
of Poisson's equation in Eq. ©. It has been verified 
that this gives a very good description of the potential 
in the realistic junction in Fig. [TJ£ The spatial profiles 
of the source-drain and gate voltages follows from the 
a-functions which are solutions to Laplace equation in 
Eq. ([5]). For simplicity, we approximate these by sim- 
ple linear functions. The coupling to the gate electrode 
is set to unity, i.e. a ga tc = 1. In realistic descriptions 
of the gate potential and in experimental realizations 
of nanoscale junctions this number is often rather low 
(~ 0.1 — 0.2)&li2rf2i and may vary over the spatial ex- 
tend of the molecule.— However, under the assumption 
that the gate potential is constant on the molecule, the 
gate-coupling parameter only serves as a scaling factor 
for the shift of the molecular energy levels. The source- 
drain voltage is modeled by a linear ramp between the 
two electrodes with the voltage applied symmetrically to 
the left and right electrode, 

V{z) = ^-V sA ^ (8) 

where z — corresponds to the position of the left elec- 
trode and L is the electrode spacing. With the chemical 
potentials of the leads given by 

/i L = E f - eV sd /2 and fi R = E f + eV ad /2, (9) 

where Ef is the equilibrium Fermi level, the positive cur- 
rent direction is from left to right. Within a microscopic 
description of both the leads and the molecule, the align- 
ment between the lead Fermi levels and the molecular lev- 
els follows directly. However, for the simplified descrip- 
tion of the leads adopted here, the alignment is treated as 
a parameter, and the equilibrium Fermi levels are chosen 
to reside in the middle of the gap of the benzene molecule. 
The level alignment is illustrated schematically in Fig. [3Jd 
for a situation with voltages applied to the electrodes. 



B. Validity of an electrostatic description 

The classical electrostatic treatment of the junction en- 
vironment in the Hamiltonian in Eq. (TTJ) relies on certain 
assumptions. First of all, the dimensions of the electrodes 
and gate dielectric need to be larger than the screening 
lengths in the respective materials in order for a classi- 
cal description to be valid. Secondly, the classical treat- 
ment of the environment assumes that there is no sig- 
nificant overlap between the quantum mechanical region 
and the environment. In the Coulomb blockade regime, 
this seems to be a fair assumption. Last, since the en- 
vironment is described with electrostatics, the time scale 
of the electrons on the molecule needs to be slower than 
the response times of the metallic electrodes and gate 
dielectric. This requires the hopping integrals t on the 
molecule to be smaller than typical values for the plas- 
mon and phonon energies. With molecular level spacing 
being on the order of eV for conjugated molecules, this is 
not necessarily the case. Nevertheless, the current carry- 
ing electrons have residence times on the molecule much 
longer than the response time of the environment. So the 
question is not whether or not the image charge effect 
should be included, but rather if it should be treated as 
an instantaneous interaction as in Eq. (J3J or at the level 
of mean- field theory where the environment only sees the 
mean occupation of the electrons on the molecule. In the 
following we pursue the first direction. We have verified 
that a mean-field treatment at the level of Hartree-Fock 
does not change the main conclusions of the present work. 

C. Current 

We shall here focus on the weak-coupling regime where 
the occupation probabilities of the molecular states and 
the current can be obtained within a master equation ap- 
proach. In the case of degeneracies between the molec- 
ular states, a master equation for the density matrix 
that retains the coherence between the degenerate states 
must be considered] 13 > 21 However, it turns out that in the 
present study such degeneracies are broken by the inter- 
action with the image charges in the junction environ- 
ment, which leaves the simpler rate-equation approach 
valid, so that only the diagonal elements of the reduced 
density matrix need to be retained i 25 i 26 

To lowest order in the coupling between the molecule 
and the leads in Eq. @, the transition rate between the 
i'th TV-electron state and the j'th N + 1-electron state 
due to tunneling from lead a is given by Fermi's golden 
rule 

HT +llJ = ^■i^{N)f a {E ij ), (io) 

where f a is the Fermi distribution for the lead electrons, 
Eij = E^ +l — Ef corresponding to the molecular ion- 
ization energies and electron affinities, and 

l"UN) = Pa t 2 a \ (N + l,j\cUN,i}\ 2 (11) 
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is the product of the lead density of states p a , the tunnel 
coupling t a , and the transition matrix element between 
the two molecular states involved in the addition of an 
electron to the molecule from lead a. For the opposite 
process where an electron tunnels from the molecule to 
the drain electrode, the transition rate is given by 



2-7T 

_ 7 -(jV-l)( l-f a (E„ 



(12) 



where the fermi factor ensures that there is an empty 
state in the lead. We here restrict our discussion to the 
case of identical left and right leads, modeled by normal 
metals with a constant density of states. The presence of 
a fast (on the scale of the time between tunneling events) 
energy relaxation mechanism in the leads justifies the 
equilibrium description of the lead electrons in Eqs. (fTO"]) 
and (fT2]) . 

The rate-equations for the occupations of the molecu- 
lar states now reads 



P N ,i = £ 



(13) 



Together with the normalization condition ^2 N i Pn.i = 1 
the rate-equations can be solved in steady-state, i.e. 
Pna — 0, for the occupation probabilities Pn.i- From 
the steady-state occupations, the current through the 
molecule follows by evaluating the total rate of electrons 
from lead a 



N,ij 



-fW.ifsi+i.j 



PN.i^ N— 1,3 



(14) 



for a — L/R, respectively. 

The main objective is here to demonstrate how the 
inclusion of the image charge effect affects the molecular 
states and thereby also the IV characteristic of the molec- 
ular junction. For this purpose, the energies and states 
appearing in Eqs. (|TU)) . (fTTj) . and (TT21 must be calcu- 
lated on the basis of the junction Hamiltonian in Eq. ([3)) 
which fully accounts for the interaction with the image 
charges. The extend to which specific molecular prop- 
erties are affected is highly dependent on the molecule 
of interest and its configuration in the nanojunction. A 
general discussion is therefore not possible, why the rest 
of the paper seeks to address the importance of the image 
charge effect in the often studied benzene SET i 11 ' 20 ^ 1 



is set to 
A1 2 3 . 



10 corresponding to the high-K dielectric 



In experimental situations both the position of the 
molecule and the electrode spacing are highly uncontrol- 
lable parameters. We therefore consider different setups 
where the position and orientation of the molecule with 
respect to the electrodes are varied. In the symmetric 
setup shown in Fig. [Tp, the molecule is placed in the 
middle of the junction with the two end atoms facing the 
electrodes^ In this reference setup, the distance between 
the hydrogen atoms facing the electrodes at the ends of 
the molecule and the electrostatic boundaries of the elec- 
trodes is set to 1.2 A. With the so-called image charge 
plane lying ~ 1.0 A outside the outermost atomic layer of 
a surface, — this corresponds to a weak bond between the 
surface atoms of the electrodes and the benzene molecule 
(see also App. |A"|) . To account for the longer range of the 
p z -orbitals on the carbon atoms, the distance to the gate 
dielectric is chosen slightly larger to 2 A. In the following, 
all changes will be with respect to the above described 
reference setup. 

In order to model experimentally more relevant situa- 
tions, we also consider the following two setups and com- 
binations hereof. In the first, the molecule is rotated by 
an angle 6 with respect to its symmetric setup in Fig. [Tp 
around its six-fold rotational symmetry axis perpendic- 
ular to the plane of the molecule. In the second, an 
asymmetric setup where the distance to one of the elec- 
trodes is increased to twice the distance in the reference 
setup is considered. This leads to a smaller image charge 
interaction with the most distant electrode. While the 
symmetric setup corresponds to electrode couplings in 
the para configuration (marked by red dots in Fig. [TJ, 
coupling also in the meta configuration (green dots in 
Fig. [TJ is likely to occur in the asymmetric setup where 
the relative difference in the distance between the dis- 
tant electrode and the two coupling sites is smaller. The 
combination of the above described situations is shown in 
Fig. [TH where the molecule is placed in a rotated position 
closer to the left electrode. 



A. Pariser-Parr-Pople Hamiltonian 



III. BENZENE SET 

In the remaining of the paper we consider the transport 
through the benzene SET illustrated in Fig.[JJi. It has the 
molecule lying flat on a gate oxide between the source and 
drain electrodes. The dielectric constant of the gate oxide 



In conjugated molecules the sp 2 hybridization of the 
carbon atom results in an energy separation between the 
bonding c-orbitals and the higher lying 7r-orbitals which 
have mainly p z character. Quantitative predictions of the 
low energy excitations can therefore be obtained with the 
simple Pariser-Parr-Pople 2 ^— (PPP) description. The 
PPP Hamiltonian which includes only the 7r-electron sys- 
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tern is given by 



(15) 



where = rij-f + and Z.; = 1 are the occupation 
and valence, respectively, of the p 2 -orbital at the car- 
bon site i of the molecule. Apart from hopping between 
the p z orbitals on neighboring sites in the second term, 
the PPP description also includes onsite and long-ranged 
Coulomb interactions between the 7r-electrons. Due to 
the effective 7r-electron description, the parameters U and 
V cannot be identified with the usual matrix elements 
of the Coulomb interaction since screening effects from 
the CT-electrons reduce the Coulomb interaction. Here, 
the Ohno parametrization^ 2 - is used for the long-ranged 
Coulomb interactions, 



U 



(16) 



where a — [7/(14.397 eV A) and is the distance be- 
tween the carbon sites i and j. It ensures that the bare 
Coulomb interaction is recovered for large distances, i.e. 
Vij — > 1/r as Tij — >■ oo, while for shorter distances the 
(screened) onsite Coulomb interaction U is approached. 
The remaining parameters of the PPP Hamiltonian in 
Eq. ([T5|l have been chosen as £j = eV, t — 2.539 eV 
and U = 10.06 eV. This set of parameters have been fit- 
ted to experimental excitation energies and hence gives a 
quantitative description of the excitations in the benzene 
molecule^ We note that in previous studies of Coulomb 
blockade transport through the benzene molecule ; 11 ' 2 ^ 21 
different sets of parameters have been used. The ex- 
citation spectrum of the isolated benzene molecule re- 
ported here, might therefore differ slightly from these 
other works. It should be stressed that despite the large 
value of the onsite Coulomb interaction U, the inclu- 
sion of long-ranged Coulomb interactions Vij in the PPP 
description results in weakly correlated states that to a 
good approximation can be described with Hartree-Fock 
theory <2i 



The PPP Hamiltonian in Eq. (115)) describes only the 
isolated molecule. In order to include the interaction with 
the junction environment, the terms in the Hamiltonian 
in Eq. ^ involving the G part of the Green's function 
and the applied voltages Vi must be taken care of sepa- 
rately. The conversion of these terms from the real-space 
representation in Eq. to the atomic p z -orbital basis 
of the PPP Hamiltonian in Eq. (TT5|) is given in App. [B] 
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FIG. 2: (Color online) Matrix elements (in eV) of the 
Coulomb interaction, given by the Ohno parametrization 



in Eq. (|16p (left), and the image charge interaction Vij in 
Eq. (|17[) (right). The numbers on the axes denote the site 
indices of the carbon atoms in the benzene molecule with site 
1 and 4 corresponding to the carbon sites marked with red 
dots in Fig. [lb. The large value of the matrix elements of the 
image charge interaction results in a significant renormaliza- 
tion of the benzene Hamiltonian and a breaking of its high 
six-fold rotational symmetry. 



This leads to the following additional terms 



^mol-cnv F H cnv ^ 



vr + \vu + v? 



(17) 



and Vij are the matrix elements of the 



where V/ on , V^ ext 
induced junction potential Vi on due to the ionic cores, the 
matrix element of the external potential V ex t = Si a iVi 
from the applied voltages Vi , and the two-particle matrix 
elements of G in the basis of the atomic p z -orbitals, re- 
spectively. In order to adapt to the level of complexity of 
the PPP Hamiltonian, only the diagonal matrix elements 
of the terms in Eq. (fT7|) are retained here. This should 
not affect the main conclusions of the paper. 

As is evident from Eq. p7)l . the interaction with the 
image charges leads to a renormalization of both the on- 
site energies and the Coulomb interaction. The two first 
terms inside the square brackets stem from the image 
charge of the ions and the image charge of the electron 
itself. The third term inside the brackets is the shift 
in the onsite energies due to the applied voltages. The 
two last terms in the last line account for the interac- 
tion with the image charges from all the other electrons. 
These terms have the usual form of the Coulomb inter- 
action and thus correspond to a renormalization of the 
interactions of the bare molecule in Eq. (fT6l) . 

The matrix elements of the Coulomb ^interaction in 
Eq. (jT5J) and the image charge interaction Vij (see App.lBl 
for details) are shown in Fig. [2] for the symmetric setup. 
As can be seen, the matrix elements of the image 
charge interaction is on the order of several electron- 
volts and therefore result in a strong renormalization of 
the Coulomb interactions on the molecule. With the 



large renormalization of both the onsite energies and 
the Coulomb interaction in the molecular Hamiltonian 
in Eq. (|15|) . its is not surprising that the image charge 
effect has a considerable impact on the molecular states, 
their energies and symmetry. 

In the following the many-body Hamiltonian given by 
the sum of the contributions in Eqs. (fT5| and (|17p is diag- 
onalized directly in the Fock space of many-body states. 
Since the Hamiltonian commutes with the number oper- 
ator N — J2ia ™io an d the ^-projection S z of the total 
spin, each of the (N, S z )-subblocks of the Fock space are 
diagonalized separately. For the neutral N = 6 state of 
the molecule, the dimension of this subblock is 400 x 400 
implying that the Hamiltonian can be diagonalized with 
standard diagonalization routines. This yields the many- 
body states 

|7V,z)=5>|^) (18) 

n 

and their corresponding energies , where {\<fin)} de- 
notes the possible iV-electron configurations with spin 
S z and c n are the expansion coefficients. For example, in 
the singlet ground-state of the neutral isolated molecule 
where the Hamiltonian is given by Eq. (|15|) . the config- 
uration \(t>n =6 } = with the spins aligned op- 
positely on neighboring sites is the one with the largest 
weight. 

In the following, we will focus on transport through 
the benzene molecule at positive gate voltages. This 
corresponds to the situation illustrated in Fig. [3Jd where 
the affinity level is located in the bias window. Due to 
the electron-hole symmetry of the PPP Hamiltonian, the 
transport through the positively charged cation will be 
identical. The states relevant for the transport at posi- 
tive gate voltages are illustrated schematically in Fig.[5h.. 
This includes the ground-state of the neutral molecule 
(N = 6) and the degenerate ground-state of the singly- 
charged anion (N = 7) which is split up by the image 
charge effect as illustrated in the right plot. The next 
excited states lie more than 2 eV above these states and 
will hence not be active under moderate source-drain bi- 
ases. 



B. Addition energy and breaking of degeneracies 

The most apparent effect of the interaction with the 
image charges is a strong renormalization of the molec- 
ular charging energies which has been observed both 
experimentally 6 -^ and theoretically £2 This results in a 
large reduction of the addition energy, given by the dif- 
ference between the ionization potential (IP) and electron 
affinity (EA), 

£ add = IP - EA = < +1 + Eg' 1 - 2Eg, (19) 

compared to its value for the isolated molecule. As shown 
in the first column of Tab. HI the addition energy of the 
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a) Isolated SET 




FIG. 3: (Color online) Schematic overview of the states and 
transport characteristics of the benzene SET. a) States of 
the benzene molecule important for the low-bias transport 
through the molecule. In the SET environment (right), the 
image charge effect lifts the degeneracy of the anionic ground- 
state of the isolated molecule (left) producing a splitting A 
between the two states labeled \s) and \a) (see text). The 
transition rates F between the states are indicated with ar- 
rows. The transition marked with the (red) dashed arrow is 
responsible for the occurrence of NDR. b) Alignment between 
the molecular levels and the the Fermi levels of the leads in 
an out-of-equilibrium situation, c) Current as a function of 
source-drain voltage. With state \a) being a so-called blocking 
state with a slow exit rate illustrated by the dashed red arrow 
in a) , it is likely to produce an NDR feature when introduced 
in the bias window. 



benzene molecule is reduced by up to ~ 3.6 eV in setups 
similar to those in Fig. [TJ This is of the same size as the 
reductions found with ab-initio descriptions of the ben- 
zene molecule in similar environments ! 35 ' 36 Due to the 
smaller image charge effect from the most distant elec- 
trode in the asymmetric setup, the reduction of the addi- 
tion energy is smaller in this case. Similar reductions of 
the molecular HOMO-LUMO gaps are expected to occur 
in single-molecule junctions with a stronger coupling to 
the leads where coherent transport dominated 

The large renormalization of the molecular energy lev- 
els is also manifested in the charge distribution on the 
molecule where the attractive nature of the oppositely 
charged image charges polarizes the molecule. The mag- 
nitude of this effect can be quantified in terms of the site 
occupations on the benzene ring which follow from the 
expectation value (hi) — (N\fii\N). In the symmetric 
setup in Fig.[TjD, we find that ~ 0.3 of the added electron 
in the anion resides on each of the end atoms closest to 
the electrodes while only ^0.1 resides on each of the four 
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E a dd 


A 


Symmetry 


Isolated 


11.38 




D 6h 


SET (symmetric) 


7.77 


72 


D 2h 


SET (asymmetric) 


8.01 


43 


C2v 


SET (rotated) 


7.73 


76 





TABLE I: Addition energy (in eV) and splitting A (in meV) 
of the degenerate ground-state of the negatively charged anion 
(N = 7) in different situations. The symmetric and antisym- 
metric configurations correspond to the setups illustrated in 
Figs.QJ) and c, respectively. In the rotated setup, the molecule 
has been rotated by an angle 6 = tt/6 around its six- fold ro- 
tational symmetry axis. The reduction of ~ 3.6 eV for the 
addition energy in the symmetric and rotated setups is of the 
same size as the reductions found with ab-initio GW/DFT 
calculations for a benzene molecule on a graphite surface^ 
and in a SET environment.— The smaller reduction observed 
for the addition energy in the third row is a consequence of the 
asymmetric setup which reduces image charge effect from the 
more distant electrode. In the last column the point group 
symmetry of the Hamiltonian in the absence of an applied 
bias voltage is given. 



center atoms. For the exact diagonalization of the Hamil- 
tonian given here, this charge rearrangement corresponds 
to a change in the weights c n of the different iV-particle 
configurations in the many-body state compared to the 
isolated molecule. 

In transport measurements, the addition energy can 
be inferred from the height of the Coulomb diamonds 
in the charge stability diagram. 38 However, even when 
taking into account the image charge effect, the addi- 
tion energy of the benzene molecule is large compared 
to experimentally accessible source-drain voltages, and 
hence the observation of full Coulomb diamonds for such 
a small molecule seems unlikely. We will therefore fo- 
cus on the IV characteristics at lower biases, where the 
image charge effect leaves its fingerprint in the form of 
an additional molecular level that results from a broken 
symmetry in the molecule. 

The PPP Hamiltonian for the isolated benzene 
molecule in Eq. (1X3]) belongs to the Dgh point group. The 
symmetry and degeneracies of the different charge states 
of the PPP Hamiltonian for benzene have been consid- 
ered in detail in Refllll. The symmetry of the ground- 
states for the neutral molecule and the anion is Ai g and 
, respectively, with the latter having a twofold orbital 
degeneracy on top of its spin degeneracy. The states are 
illustrated schematically in Fig. [3ji. 

When taking into account the image charge effect, 
the symmetry of the full junction Hamiltonian given by 
Eqs. (fl"5|) and (TP?)) is reduced with respect to that of the 
of the isolated molecule. In this case, the symmetry of the 
junction Hamiltonian reflects the symmetry of the com- 
bined molecule plus junction setup. The point groups 
of the Hamiltonian for the different setups are listed in 
the last column of Tab. HI For example, in the sym- 



metric setup in Fig. [TJa, the symmetry is reduced to the 
Dih point group. As a consequence, the orbital degen- 
eracy of the anion E 2u ground-state is lifted resulting in 
a splitting A between the symmetry-broken states. The 
situation is illustrated schematically in the right part of 
Fig. (the labeling of the states is explained in the next 
section). The splitting of the degenerate anion state is 
listed in the second column of Tab. |T] for the different 
setups. In all cases, the splitting is of considerable size 
(40 — 80 meV). In the regime k^T > T considered here, 
the splitting thus exceeds the level broadening Y even at 
room-temperature in the Coulomb blockaded junctions. 
With the distance to the image plane of the electrodes 
in the symmetric setup increased to ~ 4 A, the splitting 
remains on the order of 10 meV. Hence, irrespective of 
the exact alignment between the molecule and the elec- 
trodes, the splitting is large enough for the split states 
to appear as individual resonances in the charge stability 
diagram at sufficiently low temperatures. 

In the junction considered here, the electrodes and the 
gate dielectric affect the molecular symmetry differently 
Since the molecule is lying flat on the gate dielectric, 
it does not break the symmetry of molecule. The image 
charge effect from the electrodes is therefore most impor- 
tant for the observed lifting of the degeneracies. On the 
other hand, both the gate dielectric and the electrodes 
contributes equally to the reduction of the addition en- 
ergy 3 



C. Selection rules 

The charge transport through the molecule is to a high 
degree determined by the transition matrix elements in 
Eq. pip . For molecules with symmetries, group theoreti- 
cal arguments can be used to derive selection rules for the 
transition matrix elements between the involved states. 
The selection rules for the isolated benzene molecule have 
been considered in detail in RefflTI. Due to the lower sym- 
metry of the full junction Hamiltonian considered here, 
the following analysis differs slightly. 

For the symmetric setup in Fig.[T|3 the Hamiltonian be- 
longs to the D2/1 point group. The elements of the point 
group and the Hamiltonian therefore posses a common 
set of cigenstates. Here, the symmetry of interest is the 
symmetry operation a v which is a reflection about the 
plane through the two carbon atoms closest to the elec- 
trodes and perpendicular to the molecular plane. The 
eigenstates can therefore be classified by the eigenvalues 
±1 corresponding to symmetric and antisymmetric states 
with respect to reflections in the mirror plane a v . The 
ground-state |0) of the neutral molecule is symmetric, 
i.e. a v \0) = |0). For the split ground-state of the an- 
ion, the lower and higher lying states are symmetric and 
antisymmetric, respectively. We therefore label them as 
in Fig. [3] where \s) denotes the symmetric and \a) the 
antisymmetric excited state of the anion. 

We now consider coupling in the para configuration 
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FIG. 4: (Color online) Transition matrix elements 
I (N + 1, s/a|cJ, CT |iV, 0) I between the ground-state |0) of the 
neutral molecule (N = 6) and the symmetric/antisymmetric 
state |s)/|a) of the anion (N = 7) vs angle of rotation of the 
benzene molecule. The rotation angle is with respect to the 
symmetric setup in Fig. [lb. The red and green lines refer to 
coupling at the para and meta sites denoted by red and green 
dots in Fig. [TJ respectively. 



and meta (green lines) site are shown. As evident from 
Fig. [TJ the para coupling sites become the meta site and 
vice versa under a rotation of 9 = tt/3. For 9 = tt/6 
the two coupling sites are equivalent. This is reflected in 
the mirror symmetry between the red and green lines in 
Fig.@]which meet at 9 — tt/Q. 

As discussed above, the para matrix element vanishes 
for the antisymmetric state at 9 = 0. The other tran- 
sition matrix elements all have finite values in the non- 
rotated setup. For 9 =/= the mirror symmetry of the 
junction is broken and the selection rules derived from 
Eq. (|2T)1) no longer apply. As a result, the transition ma- 
trix element for the antisymmetric state with coupling at 
the para site acquires a finite value. Under these circum- 
stances, transport via the otherwise dark state becomes 
possible. As will be discussed in further detail below, 
the small magnitude of the associated transition matrix 
element for the antisymmetric state in Fig. [T] is likely 
to cause NDR. The occurrence of NDR from a low-lying 
symmetry-split state may therefore serve as a fingerprint 
of a broken symmetry in the molecule. 



D. Current and NDR 



where the leads couple to the two atoms facing the elec- 
trodes in Fig. [TJ). Since the coupling atoms lie in the 
mirror plane a v , the transition matrix element between 
the symmetric N = 6 ground-state and the antisymmet- 
ric N = 7 state can be shown to fulfill the following 
equality 

(A+l,a|cL|A,0) - -(N + l,a\ci a \N,0), (20) 

implying that the matrix element vanishes^ Transitions 
to the antisymmetric state are therefore forbidden which 
results in a so-called dark state that cannot be observed 
in transport measurements. For transitions between the 
two symmetric states no such restriction exist. In other 
words, this is a statement that the operator c^ CT preserves 
the symmetry of the state when coupling in para configu- 
ration. In meta configuration, where one of the couplings 
is shifted to the neighboring site, the creation operator 
for the shifted coupling site no longer preserves the mirror 
symmetry of the states. Hence, the matrix element will 
be non-zero for both the symmetric and antisymmetric 
state of the anion. 

Since the mirror symmetry a v is an element of the C<i v 
point group, the selection rules derived above also apply 
in the asymmetric setup in Fig. [TJ;. One way to break 
the mirror symmetry, is to rotate the benzene molecule 
around its six-fold rotational symmetry axis as sketched 
in Fig.QJl. This breaks all the symmetries in the molecule 
and destroys the selection rules for the transition matrix 
elements. Figure [4] shows the absolute value of the tran- 
sition matrix elements as a function of the rotation angle 
9 with respect to the symmetric setup in Fig. [TJ5. The 
matrix elements for coupling to both the para (red lines) 



At relatively low biases, only the low-lying states illus- 
trated in Fig. [3Ji will be active in the charge transport 
through the molecule. More specifically, only the ground- 
state of the neutral molecule and the split degenerate 
ground-state of the singly charged molecule needs to be 
considered. Given that the temperature is low enough 
and k-gT > T, the large splittings A of the anion states 
in Tab.UJbring the junction outside the quasi-degenerate 
regime T ~ A where coherence between the degenerate 
states is important! 11 ! 13 Hence, the out-of-equilibrium oc- 
cupations of the molecular many-body states and the cur- 
rent can be obtain with the conventional rate-equation 
approach described in Sec. Ill CI 

With only three states — two of which are the 
symmetry-split doublets of the singly charged molecule — 
participating in the transport the stationary rate- 
equations take the simple form 





~(r s o + r a o) 
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{Pa) 



= (21) 



where the 0, s and a subscripts denote the ground-state of 
the neutral molecule and the symmetric ground-state and 
first excited antisymmetric state of the charged molecule, 
respectively. The rates = Tfj +Tfj have contributions 
from tunneling to both the left and right leads. For the 
states of the charged molecule P 8 / a denotes the occupa- 
tions of the individual spin up and down states of the 
doublet which are equal, i.e. P^ = = Pi. Together 
with the normalization condition Pq + 2P S + 2P a = 1 the 
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rate-equations have the solution 
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(22) 



where T 2 = r 05 .r 0a + 2(r a0 r 0s + r s0 r 0a ). The current 
through the molecule, here evaluated at the left lead, is 
then given by 

/ = -2e [p (rf + r^ ) - p s r^ s - p a r% a ] 

j^R t j^L i t^R j^L -nR -n j^L j^R -n j^L 
1 sO 1 Oai Os T 1 aO 1 Os 1 Oa — 1 Os 1 Oa^ S Q ~ 1 Oa 1 Os^ a O 



= 2e 



r 2 



(23) 



where the factor of 2 comes from the spin degeneracy of 
the anionic doublet states. 

In the following we wish to establish a general condi- 
tion for 7-factors in Eq. ([TT|) under which NDR occurs 
for the three state system of the symmetry-broken ben- 
zene molecule. Similar considerations have been given in 
for a spinless three state system. We here gener- 
alize this condition to take into account the spin degen- 
eracy of the two anionic states. The result is completely 
general and can be applied to similar three state systems 
which generally occur in single-molecule junctions when 
an excited state of a charged molecule becomes accessi- 
ble^ 

In the case of NDR, the current decreases at the voltage 
where the excited state enters the bias window. The 
situation is illustrated schematically in Fig. [3fc. Here, 
I s is the current when only the level for the symmetric 
state of the anion is located in the bias window. For 
larger values of the bias voltage where the antisymmetric 
state becomes accessible, the current either increases or 
decreases to the value I sa . In this case, the current is 
given by the expression in Eq. (|23j) . In the former case, 
the Fermi factor in Eq. (|10p is zero for the antisymmetric 
state and the expression for the current simplifies to 



/ = -2e [p rf - p s r L s ] 

yR ~pL ~pR pL 

_ 2g _ sO Os 1 0s i s0 



Os 



2r s o 



(24) 



The NDR, illustrated by the (red) dashed line in Fig. [3fc, 
occurs when \I S \ > \I sa \- At the horizontal plateaus in 
between the jumps in the current, the distance \E s / a —[j, a \ 
between the molecular levels and the chemical poten- 
tials of the leads is assumed to be larger than the ther- 
mal energy k^T. The Fermi factors f a and 1 — f a in 
Eqs. (flO|) and (fT2j) for the rates can therefore be set to 
either unity or zero. Together with the current expres- 
sions in Eqs. (|23l) and ([2"4")l . the inequality for the currents 
above leads to the following conditions for the 7-factors 
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(25) 



which must be fulfilled in order to have NDR at positive 
or negative bias voltage, respectively. These inequalities 
can be fulfilled when the excited antisymmetric state of 
the anion is a so-called blocking state which has a small 
exit rate, i.e. 7g a must be small for the drain electrode 
compared to at least one of the 7-factors for the symmet- 
ric state. The existence of such a small exit rate results 
in a slowing down of the charge transfer dynamics when 
the corresponding state enters the bias window and is 
therefore often accompanied by NDR. In Fig. the cor- 
responding transition is marked by the (red) dashed ar- 
row. It should be noted that the inequalities above only 
require the exit rate for the excited state to be small and 
put no constraint on the in rate. 

A situation with a small exit rate can be realized in 
either of the following two ways (or both). In the case of 
strongly asymmetric tunnel couplings, e.g. <C tn, the 
exit rate becomes small when the coupling to the drain 
electrode is the weaker. However, given that the tran- 
sition matrix elements for the states in Eq. (TTTj) are the 
same, this situation results in a small exit rate for both 
the symmetric and antisymmetric state and is therefore 
not sufficient for the occurrence of NDR. A difference 
in the transition matrix elements between the states in 
Eq. (|lip on top of the asymmetry in the tunnel cou- 
plings is therefore required to fulfill the inequalities in 
Eq. (1251) . On the other hand, for symmetric tunnel cou- 
plings, i.e. — tR, NDR is also possible. Here, the 
inequalities Eq. (|25[) may be satisfied by the transition 
matrix elements alone when the matrix element for the 
excited state is small. A recent theoretical study has 
demonstrated that this property can be designed into 
conjugated molecules by functionalizing them with dif- 
ferent types of chemical groups 4°- This results in asym- 
metric molecular orbitals for the excited states and con- 
sequently a small transition matrix element at one end 
of the molecule that produces NDR. As discussed in the 
previous section, the small transition matrix element is 
here provided by the destroyed selection rules in the low- 
symmetry setup in Fig. [T}i where the molecule is placed 
in a rotated configuration. 



IV. IV CHARACTERISTICS AND STABILITY 
DIAGRAMS 



In the following section the IV characteristics of the dif- 
ferent setups is considered. In all cases the temperature 
of the lead electrons is set to k^T = 5 meV and the equi- 
librium Fermi level of the leads is positioned in the mid- 
dle of the molecular gap. Furthermore, since only trans- 
port via the negatively charged anion of the molecule is 
considered, the discussion is restricted to positive gate 
voltages. It should be emphasized that the many-body 
states, their energy, and the transition matrix elements 
in Eq. (TTT1) are calculated for each value of the applied 
voltages. The IV characteristics reported in the follow- 
ing therefore include possible shifts of the molecular lev- 
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FIG. 5: (Color online) IV characteristic for the symmetric setup and coupling in the para configuration. Left: Stability diagram 
for 6 — n/12 ~ 0.26, k^T — 5 meV and ii = til- Right: Absolute value of the current vs gate voltage at V s d = 0.2 V for 
different values of 6. As evident from the blue 9 = curve, the selection rules for the transition matrix elements in Eq. (|11[1 
completely expel the antisymmetric anionic state from participating in the transport and no NDR is observed. The varying 
position of the current onset and the current dip is caused by a weak ^-dependence of the molecular energy levels. 



els with the applied bias voltage. In particular, different 
capacitive couplings to the source and drain electrodes 
arising from the asymmetric setup will result in IV char- 
acteristics that are asymmetric in the bias voltage. 



A. Symmetric setup with coupling in the para 
configuration 

Figure [5] summarizes the IV characteristics of the sym- 
metric setup in Fig. QJ>. Due to the selection rules, the 
antisymmetric state is not active in the transport for the 
non-rotated case when coupled in the para configuration. 
However, when the molecule is rotated the matrix ele- 
ment for the antisymmetric state acquires a finite value 
(see Fig. |4|) allowing the antisymmetric state to be pop- 
ulated when it enters the bias window. The left plot in 
Fig. [5] shows the stability diagram for a rotation angle 
8 = 7r/12 of the benzene molecule and symmetric tunnel 
couplings tL = tn. At this value for 8, the ratio be- 
tween the transition matrix elements for the symmetric 
and antisymmetric states with coupling to the para site 
is |(s|cJ [Cr |0)/(a|cj lo .|0)| ~ 3.7. Reverting to the condi- 
tion for the occurrence of NDR in Eq. (|25p , this value for 
the ratio between the transition matrix elements is found 
to fulfill the inequalities in the case of symmetric tunnel 
couplings. Consequently, an NDR feature appears in the 
stability diagram at the voltages where the antisymmet- 
ric state becomes accessible. The right plot in Fig. [5] 
shows the absolute value of the current as a function of 
the gate voltage at V s d = 0.2 V for different values of the 
rotation angle of the benzene molecule. Due to the sym- 
metry forbidden population of the antisymmetric state 
at 8 = 0, this state remains unpopulated at the gate 
voltage where it enters the bias window. Therefore, the 
corresponding curve in the right plot of Fig. [3] does not 



show any change in the current at the corresponding gate 
voltage. For increasing rotation angles, both the NDR ef- 
fect and the current level prior to the onset of the NDR 
is seen to decrease. This trend is a direct consequence of 
the variation of the transition matrix elements shown in 
Fig. d] when the molecule is rotated with respect to the 
electrodes. Likewise, the shift of the onset for the cur- 
rent and NDR originates from small changes in the level 
positions when the molecule is rotated. 

The linear character of the Coulomb diamond bound- 
aries reveals that the bias voltage does not have any ap- 
parent effect on the molecular states and their energies. 
Only at much higher voltages may this become impor- 
tant. For longer molecules, however, where the voltage 
drop over the molecule becomes significant, a larger im- 
pact on the molecular states is possible^ These consid- 
eration also indicate that the symmetry-breaking effect 
of the bias voltage is here negligible compared to that of 
the image charge effect. 



B. Asymmetric setup with coupling in both para 
and meta configuration 

Next, we focus on the asymmetric setups illustrated in 
Fig. Q}j. Due to the larger distance to the right electrode 
it becomes relevant to address the effect of an asymmetry 
in the tunnel couplings, i.e. > tR, and coupling to the 
meta site at the right electrode. Since the values of the 
transition matrix elements in Fig.Uare left relatively un- 
affected by the asymmetry of the setup, changes in the IV 
characteristics as compared to the symmetric setup con- 
sidered in the previous section, can be attributed to the 
introduced tunnel coupling asymmetry and the change in 
the coupling site. 
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FIG. 6: (Color online) IV characteristics for the asymmetric setup where the distance to the right electrode is twice the distance 
to the left electrode and with coupling in the para configuration. Left: Stability diagram for 9 = 7r/12 ~ 0.26, &bT = 5 meV 
and tL = 10t.R. Right: Absolute value of the current vs gate voltage for different values of the tunnel coupling ratio tnjtL at 
positive and negative bias voltage. 



1. Coupling in para configuration 

Figure [6] summarizes the IV characteristics for the 
asymmetric setup with coupling in the para configura- 
tion. Here, the left plot shows the stability diagram for 
8 = ir/12 and tn/tL = 0.1. Due to the asymmetry of 
the setup, the different capacitive couplings to the left 
and right electrodes give rise to different slopes of the 
diamond edges in the stability diagram. 

A clear change in the stability diagram compared to 
the one for the symmetric setup in Fig. [5j is the disap- 
pearance of the NDR feature for positive bias voltages. 
This behavior follows from the asymmetry in the tunnel 
couplings. For positive bias voltages the electrons exit 
the molecule to the left electrode. Depending on the ratio 
between the transition matrix elements, the asymmetry 
in the tunnel couplings may result in a large exit rate, i.e. 
7l > Jr- From the NDR conditions in Eq. (f2"5")) . one sees 
that this can remove the NDR. Indeed, for tn/tL =0.1 
and the previous stated ratio of the transition matrix el- 
ements, this is the case. For these parameter values, the 
difference in the transition matrix elements for the sym- 
metric and antisymmetric states that produced NDR in 
the case of symmetric couplings, is outweighed by the 
large asymmetry in the tunnel couplings. Hence, the 
NDR feature is absent at positive bias voltages. On the 
other hand, for negative bias voltage where the NDR 
feature still appears, the exit rate for the antisymmetric 
state to the right electrode remains small and the lower 
inequality in Eq. (|25)l is fulfilled. 

Since the tunnel couplings depend exponentially on the 
distance to the electrodes, other values for their ratio 
should also be considered. The right plot in Fig. [5] shows 
current as a function of gate voltage for different val- 
ues of the ratio ia/tz, at positive (Vsd = +0.2 V) and 
negative (Vsd — —0.3 V) bias voltage. As expected, the 
current level and the dip in the current occurring at the 



position of the NDR feature in the stability diagram, are 
highly dependent on the tunnel couplings. At positive 
bias, the current dip appears only for the smallest asym- 
metry tfi/ti, — 0.5 in the couplings. From the ratio be- 
tween the transition matrix elements for the symmetric 
and antisymmetric state given in the previous section, 
the value of the ratio tii/ti, at which no dip in the cur- 
rent is observed can be deduced from Eq. ([23]) . For pos- 
itive bias, the upper inequality in Eq. (|25[) becomes an 
equality at tjtjti, ~ 0.2. At this value for the coupling ra- 
tio, the difference in the transition matrix elements has 
been compensated for by the asymmetry in the tunnel 
couplings and no change in the current occurs when the 
antisymmetric state enters the bias window. For smaller 
values of the coupling ratio, i.e. larger asymmetry in the 
tunnel couplings, the current increases slightly. Again, 
at negative bias voltage, the tunnel coupling to the right 
drain electrode is the smaller and the current dip persists 
for all the shown values of the coupling ratio. 



2. Coupling in meta configuration 

Changing the coupling site at the right electrode to 
the meta site changes things drastically. Figure [7] sum- 
marizes the IV characteristics for this case. Again, the 
left plot shows the stability diagram for the same pa- 
rameter values as used in Fig. |U As shown in Fig. |4j 
the transition matrix elements for coupling to the meta 
site are equal at 9 = n/12 for the symmetric and anti- 
symmetric state. This has the immediate consequence 
that = 7oJ, implying that the lower inequality in 
Eq. ([23)) giving the condition for NDR at negative bias 
can never be fulfilled. This is reflected in the stability 
diagram where no NDR from the antisymmetric state is 
observed at negative bias. As opposed to the situation 
with coupling in the para configuration, the resonance 
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FIG. 7: (Color online) IV characteristics for the asymmetric setup with coupling in the meta configuration. Left: Stability 
diagram for the same set of parameters as in Fig. [6] Right: Absolute value of the current vs gate voltage for different values 
of the rotation angle 9 at tj, = 10tn (upper) and the tunnel coupling ratio tR/tL at 9 = n/12 (lower). Note the collapse of the 
current for 9 = in the upper plot. 



from the antisymmetric state now appears at positive 
bias and gives rise to an increase in the current instead 
of NDR. 



The upper right plot in Fig. [7] shows how the current 
in the vicinity of this resonance varies with the rotation 
angle 8. For 8 — a complete collapse of the current is 
observed when the antisymmetric state becomes accessi- 
ble. From the selection rules discussed in Sec. IIII CI and 
Fig. 21 we recall that the transition matrix element for 
coupling to the para site vanishes, while the matrix ele- 
ment for coupling to the meta site is finite. At positive 
bias this allows the antisymmetric state to become pop- 
ulated from the right electrode. However, since the exit 
process from the para site at the left electrode is forbid- 
den by symmetry, the molecule remains trapped in the 
antisymmetric state and the current is blocked. A sim- 
ilar collapse of the current in a benzene SET has been 
reported in Refl2(il where the blocking state was popu- 
lated via radiative relaxation from a higher lying excited 
state. For 6> 7^ the current collapse as well as the NDR 
disappears, and an increase in the current as the one in 
the stability diagram appears instead. 



As in the previous section, the disappearance of the 
NDR feature in the stability diagram is caused by the 
large asymmetry in the tunnel couplings. The ratio of 
the couplings at which the current dip vanishes, can 
be inferred from the upper inequality in Eq. (|25p . For 
8 = tt/12, we find that the condition for NDR is met 
for tjijt L > 0.41. The lower right plot in Fig. [7] shows 
the current as a function of gate voltage at V s d = +0.2 
for different values of the tunnel couplings. As expected, 
only the upper blue curve with tn/tr, = 0.5 shows NDR. 



V. SUMMARY AND CONCLUSIONS 

Using a theoretical framework developed for semicon- 
ductor nanostructures, we have studied the impact of the 
image charge effect on the molecular states and the trans- 
port in a benzene single-electron transistor operating in 
the Coulomb blockade regime. As demonstrated, the im- 
age charge effect renormalizes the charging energies and 
lifts the degeneracy of the twofold orbitally degenerate 
ground-state of the singly-charged anion of the benzene 
molecule. With the resulting splitting of the degener- 
ate states exceeding both the thermal energy k^T and 
the level broadening T from the tunnel coupling to the 
leads, this has important consequences for the low-bias 
IV characteristics of the benzene SET. In particular, the 
selections rules between the transport active many-body 
states of the molecule are destroyed, which gives rise to 
the appearance of a blocking state that leads to the oc- 
currence of NDR. From the derived NDR conditions for 
a generic three-state system, the coming and going of the 
NDR feature at different parameter values has been an- 
alyzed. It is demonstrated that the appearance of the 
NDR feature is very sensitive to asymmetries in the tun- 
nel couplings to the source and drain electrodes, to the 
bias polarity, and to changes in the coupling from the 
para site to the meta site. In experimental situations, ob- 
servations of the described transport characteristics may 
be an indication of a broken symmetry in the molecule. 

Altogether, we have demonstrated that image charge 
effects play a potentially important role, not only for the 
position of the molecular levels, but also for the molec- 
ular states and their degeneracies. As mentioned in the 
introduction, experimental studies^^ have already spec- 
ulated that image charge effects affect spin-excitations of 
molecules. With the exchange coupling in effective spin 
Hamiltonians being determined by hopping integrals and 
Coulomb matrix elements between the orbitals, a large 
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FIG. 8: (Color online) Simple junction geometry for which 
Poisson's equation can be solved analytically for the electro- 
static Green's function. The junction consists of two infinite 
parallel electrode plates where the lower half-space (y < 0) 
between them is filled with an oxide of dielectric constant e r . 
The dots illustrates the image charge solution to Poisson's 
equation given in Eq. (|A5fl . 

impact on spin states could be anticipated. 

With the present work, we have paved the way for 
more detailed descriptions of single-molecule SETs tak- 
ing into account image charge effects. In this respect, 
the electrostatic Green's function for the generic junc- 
tion geometry given in App. [X] and the inclusion of the 
additional image charge terms in the PPP Hamiltonian 
in App. [B] provide a good starting point for future studies 
of the image charge effect in single-molecule junctions. 
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Appendix A: Electrostatic Green's function for a 
generic junction geometry 

In this appendix we give an analytical solution of Pois- 
son's equation for the electrostatic Green's function in 



the simplified junction geometry shown in Fig. [8] De- 
spite its simple structure, the Green's function of this 
junction provides a good description of the potential in 
more realistic junctions as the one illustrated in Fig. [TJ 

When solving Poisson's equation in Eq. ^ , the screen- 
ing induced polarization charge of dielectric regions is 
most often accounted for by the spatially dependent di- 
electric constant £ r (r). In the following a different route 
will be taken by treating the polarization charge as a 
source term in Poisson's equation. For the homogeneous 
gate dielectric in Fig. [8j the polarization charge induced 
by a unit source charge in r' will be a surface charge a 
that resides on the interface between the dielectric and 
the vacuum region. The Green's function can therefore 
be obtained from the following Poisson equation 

- e V 2 G(r, r') = 5{r - r') + ^ (Al) 

e 

with Dirichlet boundary conditions, i.e. G — at the 
surface of the metallic electrodes. The surface charge 
on the right-hand side is located in the sz-plane cr(r) = 
u(x, z)S(y) and can be obtained as the normal component 
of the polarization P at the interface. Using the relations 
D = eoE + P and D = eo£ r E between the displacement 
field D and the electric field E in linear dielectrics, the 
surface charge can be related to the electric field via 

a = P • n = e (e r - 1)E ■ n, (A2) 

where the normal component (in this case the y- 
component) is evaluated immediately below the interface, 
i.e. y = CP. Expressing the electric field by the gradient 
of the Green's function, the following relation between 
the Green's function and the surface charge is obtained 

- n • VG(r, r') = 1 a(r). (A3) 

e (e r - 1) 

This relation can be used to eliminate the surface charge 
in Eq. (IA1|) . The resulting equation for the Green's func- 
tion is solved by expanding in plane-waves and sines as 

G(r,r')= / f I ^]>>n^V^G(r^,r'). 

J IT J 7T n 

(A4) 

After much algebra, the following solution for the elec- 
trostatic Green's function is found 
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(A5) 



where r, r' belong to the vacuum region of the junc- 
tion. The solution has the intuitive image charge in- 
terpretation illustrated in Fig. [51 The analytic solution 
in Eq. (IA5|) allows for a direct identification of the two 
contributions to the Green's function indicated in the 
last equality. Here, the direct Coulomb interaction is 
given by the first term inside the square brackets for 
a = t = +1, while the remaining terms give the induced 
potential G. By comparing to finite element solutions 
of Poisson's equation ([6]) in junction geometries as the 
one illustrated in Fig. [TJ we found that the Green's func- 
tion for the simplified junction considered here to a high 
degree resembles that of more realistic junctions. 

As mentioned in the main text the distances between 
the atoms of the molecule and the electrostatic bound- 
aries of the junction must be chosen with care. The 
reason for this is that the positions of the electrostatic 
boundaries between metallic/dielectric regions and the 
vacuum region where the molecule resides do not corre- 
spond to the actual positions of the atomic surfaces in 
the junction. Typically, this so-called electrostatic image 
plane of the atomic surfaces lie ~ 1 A outside the outer- 
most atomic layer.— The distance between the atoms of 
the molecule and the surface atoms are therefore larger 
than the chosen distance to the respective image planes. 



Appendix B: Image charge Hamiltonian in a 
localized basis 

In this appendix the derivation of the PPP representa- 
tion in Eq. (TTTl) of the image charge related terms in the 
junction Hamiltonian in Eq. is outlined. 

Within the PPP description, the Hamiltonian is ex- 
pressed in the basis of the p z -orbitals on the carbon 
atoms. In this basis the terms related to the image charge 
effect in Eq. (JT7J> read 

^mol-cnv ~^ Hqiyv ~ ^ ^ V{j ^ia^ja ^ ^ ^ia^ij^ja' • 

ij,cr ij : <ra' 

( B1 ) 

Here, Vl? n (> 0) and Vij (< 0) are the matrix elements 
of the image charge potential from the ionic cores and 
the two-electron integrals of the image charge interaction 
between the electrons. The two types of matrix elements 



are given by 

V$* = |dr0*(r)Vion(r)^(r) (B2) 

and 

% = JdvJ dr%(r)| 2 G(r,r')|<M r ')| 2 , (B3) 

respectively, where the induced potential from the ions 
is given by 

Vion(r) = -y"dr'G(r,r')pion(r') 

«-^G(r,R,). (B4) 

i 

In the last equality, the charge distribution of the ionic 
cores have been approximated by (5-functions located at 
the positions Ri of the carbon atoms in the molecule. No- 
tice that the image potential from the positively charged 
ions lifts the onsite energy of the electrons. The ma- 
trix elements of the image charge interaction between 
the electrons in Eq. (IB3|) are negative and therefore cor- 
respond to a screening of the direct Coulomb interactions 
between the electrons. 

Here we adopt the level of simplicity of the PPP Hamil- 
tonian and neglect off-diagonal matrix elements of the 
first term in Eq. (|B1[) and include only the direct matrix 
elements of the image charge interaction in the second 
term. This results in the Hamiltonian in Eq. (I17|) with 
the additional terms from the applied voltages treated 
similarly. 

The matrix elements in Eqs. (|B2[) and (|B3|) are evalu- 
ated by approximating the absolute square of the orbitals 
by 5- functions centered at the atomic sites, 

K on = Jdv |^(r)| 2 T4on(r)«Mon(Ri) (B5) 

and 

% = Jdv J dr' |<Mr)| 2 G(r,r')|^(r')| 2 « G(R,,Rj), 

(B6) 

respectively. The matrix elements Vij obtained for the 
junction geometry described in the main text, are shown 
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in Fig. [2] Also shown are the matrix elements of the di- of the image charge interaction leads to significant screen- 
rect Coulomb interaction given by the Ohno parametriza- ing of the interactions on the molecule, 
tion in Eq. (fl6|) . As the figure shows, the matrix elements 
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